Brain regions downstream

Setup

library(magrittr)
library(dplyr)
library(ggplot2)
library(cowplot)
library(ggrepel)
library(ggpubr)
library(scHelper)
library(biomaRt)
library(pheatmap)
library(SummarizedExperiment)
library(qs)
system("sudo apt-get update\nsudo apt-get install libnetcdf-dev -y")
library(DEP)

tt <- Sys.time()

Load data

tmp <- qread("/work/proteomics/02_data/03_objects/preprocessed_data.qs")

df.clean <- tmp$expDesign.clean
universe <- tmp$universe
dat.clean <- tmp$dat.norm.noout # No outliers, no imputation

Metadata

# Metadata
metadata <- read.delim("/work/proteomics/02_data/Metadata_Cohort_100723.csv", sep = ";", dec = ",", header = T) %>% 
  mutate(intern = gsub("K", "", .$K.no.),
         Diagnosis = gsub("Unknown", "MSA", .$Diagnosis)) # 1418 should be MSA

dat.clean@colData$subtype <- metadata$Subtype[match(dat.clean@colData$intern, metadata$intern)]

REGION-SPECIFIC PROTEINS, NO IMPUTATION

Clean data

dat.tmp <- dat.clean

Calculate DEPs

# dat.tmp <- dat.clean
# dat.tmp@elementMetadata$ID <- dat.tmp@elementMetadata$CER_ID
# dat.tmp@elementMetadata$name <- dat.tmp@elementMetadata$CER_name
all_contrasts <- test_diff(dat.tmp, type = "all")

## Warning in test_diff(dat.tmp, type = "all"): Missing values in 'dat.tmp'

## Tested contrasts: CER_vs_FC, CER_vs_MED, CER_vs_STR, CER_vs_SN, FC_vs_MED, FC_vs_STR, FC_vs_SN, MED_vs_STR, MED_vs_SN, STR_vs_SN

## Warning: Partial NA coefficients for 253 probe(s)

dep_all <- add_rejections(all_contrasts, alpha = 0.05, lfc = 2)

Number of significant proteins

comb <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

comb %>% 
  lapply(\(x) dep_all@elementMetadata %>% 
           as.data.frame() %>% 
           .[, paste(x, "significant", sep = "_")]) %>% 
  setNames(comb) %>% 
  sapply(sum) %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(line = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "", y = "No. significant proteins", fill = "Comparison") +
  scale_fill_manual(values = ggsci::pal_cosmic()(10)) +
  geom_text(aes(y = value + 8, label = value))

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, 45)) +
  ggplot2::guides(shape = "none")

## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.

## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Correlation

plot_dist(dep_all, significant = TRUE, pal = "YlOrRd")

Clustering of DEPs per sample

ht <- plot_heatmap(dep_all, type = "centered", kmeans = F, col_limit = 3, show_row_names = FALSE,
             indicate = c("condition"))

## Warning: Missing values in 'dep_all'. Using clustering_distance = 'gower'

Clustering of DEPs per comparison

plot_heatmap(dep_all, type = "contrast", kmeans = F, 
             col_limit = 5, show_row_names = FALSE)

## Warning: Missing values in 'dep_all'. Using clustering_distance = 'gower'

Volcano plots

combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)

## $V1

## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: ggrepel: 28 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V2

## Warning: Removed 238 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: ggrepel: 140 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V3

## Warning: Removed 50 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: ggrepel: 129 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V4

## Warning: Removed 76 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: ggrepel: 159 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V5

## Warning: Removed 234 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: ggrepel: 83 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V6

## Warning: Removed 38 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $V7

## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $V8

## Warning: Removed 236 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: ggrepel: 44 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V9

## Warning: Removed 235 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $V10

## Warning: Removed 68 rows containing missing values or values outside the scale range
## (`geom_point()`).

Extract results

res <- dep_all@elementMetadata@listData %>% 
  as.data.frame()

res.split <- combn(c("CER", "FC", "MED", "STR", "SN"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>% filter(!!sym(paste0(cc, "_p.adj")) <= 0.05, abs(!!sym(paste0(cc, "_diff"))) >= 1.5)), .)

GO on DEPs

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## 

## 

res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Get markers

all.markers <- res.split %>% 
  names() %>% 
  lapply(\(comp) {
    up <- strsplit(comp, "_vs_") %>% sget(1)
    down <- strsplit(comp, "_vs_") %>% sget(2)
    diff <- paste0(comp, "_diff")
    
    up.df <- res.split[[comp]] %>% 
      filter(!!sym(diff) > 0) %>% 
      mutate(area = up) %>% 
      dplyr::select(name, !!sym(diff), area)
    
    down.df <- res.split[[comp]] %>% 
      filter(!!sym(diff) < 0) %>% 
      mutate(area = down) %>% 
      dplyr::select(name, !!sym(diff), area)
    
    out <- rbind(up.df, down.df) %>% 
      dplyr::rename(l2fc = !!sym(diff)) %>% 
      mutate(l2fc = abs(l2fc))
    
    return(out)
  }) %>% 
  bind_rows()

markers.sum <- all.markers %>% 
  dplyr::group_by(area, name) %>% 
  summarize(n = n(),
            mean_l2fc = mean(l2fc),
            max_l2fc = max(l2fc)) %>%
  data.frame()

## `summarise()` has grouped output by 'area'. You can override using the
## `.groups` argument.

markers.split <- markers.sum %>%  
  arrange(-n, -mean_l2fc) %$% 
  split(., area)

Plot top markers

markers.split %>% 
  lapply(dplyr::slice, seq(10)) %>% 
  lapply(pull, name) %>%
  Map(\(genes, name) plot_single(dep = dep_all, proteins = genes, type = "centered") + guides(col = "none") + labs(title = name), genes = ., name = names(.))

## $CER

## Warning: Removed 49 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $FC

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $MED

## Warning: Removed 79 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $SN

## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_point()`).

## 
## $STR

## Warning: Removed 57 rows containing missing values or values outside the scale range
## (`geom_point()`).

Onthology on top markers

markers.enrich <- markers.split %>% 
  lapply(dplyr::slice, seq(50)) %>% 
  lapply(pull, name) %>%
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

markers.enrich %>% 
  lapply(head, 20) %>% 
  lapply(knitr::kable)

## $CER
## 
## 
## |           |ID         |Description                                                                          |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|    zScore|  pvalue|  p.adjust|    qvalue|geneID                                                                                           | Count|
## |:----------|:----------|:------------------------------------------------------------------------------------|:---------|:--------|----------:|--------------:|---------:|-------:|---------:|---------:|:------------------------------------------------------------------------------------------------|-----:|
## |GO:1903241 |GO:1903241 |U2-type prespliceosome assembly                                                      |7/49      |22/6858  |  0.3181818|      44.532468| 17.347905| 0.0e+00| 0.0000001| 0.0000001|SF3B5/SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1                                                     |     7|
## |GO:0006310 |GO:0006310 |DNA recombination                                                                    |10/49     |112/6858 |  0.0892857|      12.496356| 10.405656| 0.0e+00| 0.0000015| 0.0000014|H1-3/H1-10/H1-5/H1-4/HMGB1/POLB/APEX1/TOP2B/FUS/HDGFL2                                           |    10|
## |GO:0006259 |GO:0006259 |DNA metabolic process                                                                |16/49     |391/6858 |  0.0409207|       5.727230|  8.165206| 0.0e+00| 0.0000015| 0.0000014|H1-3/H1-10/NFIA/H1-5/H1-4/HMGB1/POLB/PDS5B/APEX1/SF3B5/TOP2B/FUS/FTO/SF3B3/HDGFL2/MPG            |    16|
## |GO:0016071 |GO:0016071 |mRNA metabolic process                                                               |16/49     |445/6858 |  0.0359551|       5.032240|  7.461381| 0.0e+00| 0.0000073| 0.0000068|ALKBH5/APEX1/SRSF2/SF3B5/FUS/FTO/SNRPD3/PRPF40A/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1 |    16|
## |GO:0045934 |GO:0045934 |negative regulation of nucleobase-containing compound metabolic process              |16/49     |474/6858 |  0.0337553|       4.724361|  7.128824| 1.0e-07| 0.0000142| 0.0000132|H1-3/H1-10/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1           |    16|
## |GO:0051253 |GO:0051253 |negative regulation of RNA metabolic process                                         |15/49     |423/6858 |  0.0354610|       4.963092|  7.137621| 1.0e-07| 0.0000158| 0.0000147|H1-3/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/FUS/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1                 |    15|
## |GO:0045892 |GO:0045892 |negative regulation of DNA-templated transcription                                   |14/49     |365/6858 |  0.0383562|       5.368297|  7.275442| 1.0e-07| 0.0000158| 0.0000147|H1-3/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1                     |    14|
## |GO:1902679 |GO:1902679 |negative regulation of RNA biosynthetic process                                      |14/49     |367/6858 |  0.0381471|       5.339042|  7.247606| 2.0e-07| 0.0000158| 0.0000147|H1-3/H1-5/H1-4/DRAP1/HMGB1/HDGF/APEX1/SRSF2/CTBP2/NRIP2/UBE2I/ILF3/CCAR1/DR1                     |    14|
## |GO:0000245 |GO:0000245 |spliceosomal complex assembly                                                        |7/49      |60/6858  |  0.1166667|      16.328571| 10.116032| 2.0e-07| 0.0000160| 0.0000149|SF3B5/SNRPD3/SF3B1/SNRPD2/SF3B3/SNRPB/SNRPA1                                                     |     7|
## |GO:0006397 |GO:0006397 |mRNA processing                                                                      |13/49     |318/6858 |  0.0408805|       5.721602|  7.313706| 2.0e-07| 0.0000160| 0.0000149|ALKBH5/SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/CPSF6/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1               |    13|
## |GO:0000377 |GO:0000377 |RNA splicing, via transesterification reactions with bulged adenosine as nucleophile |11/49     |215/6858 |  0.0511628|       7.160702|  7.785577| 2.0e-07| 0.0000160| 0.0000149|SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1                            |    11|
## |GO:0000398 |GO:0000398 |mRNA splicing, via spliceosome                                                       |11/49     |215/6858 |  0.0511628|       7.160702|  7.785577| 2.0e-07| 0.0000160| 0.0000149|SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1                            |    11|
## |GO:0000375 |GO:0000375 |RNA splicing, via transesterification reactions                                      |11/49     |219/6858 |  0.0502283|       7.029913|  7.693168| 3.0e-07| 0.0000178| 0.0000166|SRSF2/SF3B5/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1                            |    11|
## |GO:0065004 |GO:0065004 |protein-DNA complex assembly                                                         |7/49      |69/6858  |  0.1014493|      14.198758|  9.347127| 5.0e-07| 0.0000284| 0.0000263|H1-3/H1-10/H1-5/H1-4/HMGB1/GTF2A2/DR1                                                            |     7|
## |GO:0008380 |GO:0008380 |RNA splicing                                                                         |12/49     |292/6858 |  0.0410959|       5.751747|  7.039110| 6.0e-07| 0.0000358| 0.0000332|SRSF2/SF3B5/FUS/SNRPD3/PRPF40A/SF3B1/SNRPD2/ILF3/SF3B3/SNRPB/SRSF3/SNRPA1                        |    12|
## |GO:0006281 |GO:0006281 |DNA repair                                                                           |11/49     |241/6858 |  0.0456432|       6.388179|  7.223428| 7.0e-07| 0.0000379| 0.0000352|HMGB1/POLB/PDS5B/APEX1/SF3B5/TOP2B/FUS/FTO/SF3B3/HDGFL2/MPG                                      |    11|
## |GO:0071824 |GO:0071824 |protein-DNA complex organization                                                     |7/49      |78/6858  |  0.0897436|      12.560440|  8.710243| 1.1e-06| 0.0000545| 0.0000506|H1-3/H1-10/H1-5/H1-4/HMGB1/GTF2A2/DR1                                                            |     7|
## |GO:0051052 |GO:0051052 |regulation of DNA metabolic process                                                  |10/49     |209/6858 |  0.0478469|       6.696612|  7.094716| 1.6e-06| 0.0000761| 0.0000706|H1-3/H1-10/H1-5/H1-4/HMGB1/SF3B5/TOP2B/FUS/SF3B3/HDGFL2                                          |    10|
## |GO:0000018 |GO:0000018 |regulation of DNA recombination                                                      |6/49      |56/6858  |  0.1071429|      14.995627|  8.920547| 2.4e-06| 0.0001069| 0.0000992|H1-3/H1-10/H1-5/H1-4/FUS/HDGFL2                                                                  |     6|
## |GO:0030261 |GO:0030261 |chromosome condensation                                                              |4/49      |15/6858  |  0.2666667|      37.322449| 11.945956| 3.0e-06| 0.0001252| 0.0001162|H1-3/H1-10/H1-5/H1-4                                                                             |     4|
## 
## $FC
## 
## 
## |           |ID         |Description                                  |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|    zScore|   pvalue|  p.adjust|    qvalue|geneID                                                                                                        | Count|
## |:----------|:----------|:--------------------------------------------|:---------|:--------|----------:|--------------:|---------:|--------:|---------:|---------:|:-------------------------------------------------------------------------------------------------------------|-----:|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission               |18/48     |495/6858 |  0.0363636|       5.195455|  8.135136| 0.00e+00| 0.0000011| 0.0000010|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/AMPH/RIMBP2/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2/NPTX1/SYNGR1 |    18|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling         |18/48     |495/6858 |  0.0363636|       5.195455|  8.135136| 0.00e+00| 0.0000011| 0.0000010|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/AMPH/RIMBP2/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2/NPTX1/SYNGR1 |    18|
## |GO:0046928 |GO:0046928 |regulation of neurotransmitter secretion     |7/48      |62/6858  |  0.1129032|      16.131048| 10.047368| 2.00e-07| 0.0000731| 0.0000641|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B                                                                        |     7|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission |13/48     |346/6858 |  0.0375723|       5.368136|  6.999889| 4.00e-07| 0.0001002| 0.0000877|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/STX1A/SYT1/SYT5/SYN1/SV2B/SYNGR1                                |    13|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling       |13/48     |347/6858 |  0.0374640|       5.352666|  6.985707| 4.00e-07| 0.0001002| 0.0000877|EPHA4/SLC4A10/CAMKV/CAMK2A/CSPG5/HOMER1/PLCB1/STX1A/SYT1/SYT5/SYN1/SV2B/SYNGR1                                |    13|
## |GO:0051588 |GO:0051588 |regulation of neurotransmitter transport     |7/48      |75/6858  |  0.0933333|      13.335000|  9.017227| 7.00e-07| 0.0001378| 0.0001207|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B                                                                        |     7|
## |GO:0007269 |GO:0007269 |neurotransmitter secretion                   |8/48      |116/6858 |  0.0689655|       9.853448|  8.073488| 1.10e-06| 0.0001611| 0.0001411|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2                                                                   |     8|
## |GO:0099643 |GO:0099643 |signal release from synapse                  |8/48      |116/6858 |  0.0689655|       9.853448|  8.073488| 1.10e-06| 0.0001611| 0.0001411|CAMK2A/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2                                                                   |     8|
## |GO:0006836 |GO:0006836 |neurotransmitter transport                   |9/48      |161/6858 |  0.0559006|       7.986801|  7.531220| 1.30e-06| 0.0001682| 0.0001473|CAMK2A/CSPG5/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2                                                           |     9|
## |GO:2000300 |GO:2000300 |regulation of synaptic vesicle exocytosis    |5/48      |38/6858  |  0.1315789|      18.799342|  9.236718| 5.70e-06| 0.0006456| 0.0005654|CSPG5/SYT1/SYT5/SYN1/SV2B                                                                                     |     5|
## |GO:0099504 |GO:0099504 |synaptic vesicle cycle                       |9/48      |194/6858 |  0.0463918|       6.628222|  6.676037| 6.20e-06| 0.0006456| 0.0005654|CSPG5/AMPH/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2                                                             |     9|
## |GO:0006887 |GO:0006887 |exocytosis                                   |10/48     |252/6858 |  0.0396825|       5.669643|  6.340577| 7.20e-06| 0.0006811| 0.0005964|VSNL1/CSPG5/STX1A/WIPF3/SYT1/SYT5/SYN1/SV2B/SYN2/SYNGR1                                                       |    10|
## |GO:0023061 |GO:0023061 |signal release                               |10/48     |260/6858 |  0.0384615|       5.495192|  6.203588| 9.40e-06| 0.0008292| 0.0007261|VSNL1/CAMK2A/CSPG5/PLCB1/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2                                                       |    10|
## |GO:0045055 |GO:0045055 |regulated exocytosis                         |8/48      |157/6858 |  0.0509554|       7.280255|  6.682999| 1.09e-05| 0.0008777| 0.0007686|CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2/SYNGR1                                                                   |     8|
## |GO:0050808 |GO:0050808 |synapse organization                         |12/48     |394/6858 |  0.0304569|       4.351523|  5.752462| 1.15e-05| 0.0008777| 0.0007686|EPHA4/SYNPO/CAMKV/PLXNA1/HOMER1/PLXNA4/GPM6A/RAP2A/ICAM5/SYN1/SYN2/NPTX1                                      |    12|
## |GO:0099003 |GO:0099003 |vesicle-mediated transport in synapse        |9/48      |217/6858 |  0.0414747|       5.925691|  6.190054| 1.54e-05| 0.0010499| 0.0009194|CSPG5/AMPH/SLC17A7/STX1A/SYT1/SYT5/SYN1/SV2B/SYN2                                                             |     9|
## |GO:0048858 |GO:0048858 |cell projection morphogenesis                |12/48     |406/6858 |  0.0295567|       4.222906|  5.620535| 1.56e-05| 0.0010499| 0.0009194|EPHA4/CPNE5/CSPG5/PLXNA1/OLFM1/PLXNA4/GPM6A/NELL2/RAP2A/THY1/SYT1/NPTX1                                       |    12|
## |GO:1903305 |GO:1903305 |regulation of regulated secretory pathway    |6/48      |85/6858  |  0.0705882|      10.085294|  7.075743| 2.45e-05| 0.0015574| 0.0013637|CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B                                                                               |     6|
## |GO:0016079 |GO:0016079 |synaptic vesicle exocytosis                  |6/48      |87/6858  |  0.0689655|       9.853448|  6.976857| 2.80e-05| 0.0016857| 0.0014761|CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B                                                                               |     6|
## |GO:0017157 |GO:0017157 |regulation of exocytosis                     |7/48      |137/6858 |  0.0510949|       7.300183|  6.253319| 4.05e-05| 0.0023111| 0.0020237|VSNL1/CSPG5/STX1A/SYT1/SYT5/SYN1/SV2B                                                                         |     7|
## 
## $MED
## 
## 
## |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:--|:-----------|:---------|:-------|----------:|--------------:|------:|------:|--------:|------:|:------|-----:|
## 
## $SN
## 
## 
## |ID |Description |GeneRatio |BgRatio | RichFactor| FoldEnrichment| zScore| pvalue| p.adjust| qvalue|geneID | Count|
## |:--|:-----------|:---------|:-------|----------:|--------------:|------:|------:|--------:|------:|:------|-----:|
## 
## $STR
## 
## 
## |           |ID         |Description                                     |GeneRatio |BgRatio  | RichFactor| FoldEnrichment|   zScore|    pvalue|  p.adjust|    qvalue|geneID                                                                                              | Count|
## |:----------|:----------|:-----------------------------------------------|:---------|:--------|----------:|--------------:|--------:|---------:|---------:|---------:|:---------------------------------------------------------------------------------------------------|-----:|
## |GO:0050804 |GO:0050804 |modulation of chemical synaptic transmission    |14/47     |346/6858 |  0.0404624|       5.904071| 7.775855| 0.0000000| 0.0000239| 0.0000201|ACHE/CAMKV/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A                |    14|
## |GO:0099177 |GO:0099177 |regulation of trans-synaptic signaling          |14/47     |347/6858 |  0.0403458|       5.887056| 7.760662| 0.0000000| 0.0000239| 0.0000201|ACHE/CAMKV/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A                |    14|
## |GO:0007268 |GO:0007268 |chemical synaptic transmission                  |16/47     |495/6858 |  0.0323232|       4.716441| 7.130327| 0.0000001| 0.0000239| 0.0000201|ACHE/CAMKV/SLC17A7/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A/SLC1A2 |    16|
## |GO:0098916 |GO:0098916 |anterograde trans-synaptic signaling            |16/47     |495/6858 |  0.0323232|       4.716441| 7.130327| 0.0000001| 0.0000239| 0.0000201|ACHE/CAMKV/SLC17A7/ACE/SYT1/STX1A/HOMER1/SV2B/PRKAR2B/SYNGR1/PLCB1/SYP/CSPG5/SLC24A2/PPP1R9A/SLC1A2 |    16|
## |GO:0051050 |GO:0051050 |positive regulation of transport                |13/47     |461/6858 |  0.0281996|       4.114737| 5.751673| 0.0000082| 0.0019916| 0.0016713|ACHE/ACTN2/RAB27B/ATP2B1/SYT1/STX1A/SCN4B/HOMER1/RAB27A/PLCB1/THY1/SLC9A1/SLC1A2                    |    13|
## |GO:0010959 |GO:0010959 |regulation of metal ion transport               |8/47      |188/6858 |  0.0425532|       6.209144| 6.015794| 0.0000346| 0.0069868| 0.0058630|ACTN2/ATP2B1/ACE/SCN4B/HOMER1/DPP10/THY1/SLC9A1                                                     |     8|
## |GO:0007186 |GO:0007186 |G protein-coupled receptor signaling pathway    |10/47     |318/6858 |  0.0314465|       4.588519| 5.443152| 0.0000446| 0.0077314| 0.0064879|ACTN2/NECAB2/ADCY5/ACE/GNAL/HOMER1/PRKAR2B/PLCB1/SYP/DGKI                                           |    10|
## |GO:0043270 |GO:0043270 |positive regulation of monoatomic ion transport |6/47      |101/6858 |  0.0594059|       8.668211| 6.448941| 0.0000579| 0.0078021| 0.0065471|ACTN2/ATP2B1/SCN4B/HOMER1/THY1/SLC9A1                                                               |     6|
## |GO:1990138 |GO:1990138 |neuron projection extension                     |6/47      |101/6858 |  0.0594059|       8.668211| 6.448941| 0.0000579| 0.0078021| 0.0065471|OLFM1/CPNE5/SYT1/RASAL1/PLXNA4/PLXNA1                                                               |     6|
## |GO:0030001 |GO:0030001 |metal ion transport                             |11/47     |415/6858 |  0.0265060|       3.867624| 5.006252| 0.0000836| 0.0092611| 0.0077715|ACTN2/SLC17A7/ATP2B1/ACE/SCN4B/HOMER1/DPP10/PLCB1/THY1/SLC24A2/SLC9A1                               |    11|
## |GO:0045055 |GO:0045055 |regulated exocytosis                            |7/47      |157/6858 |  0.0445860|       6.505760| 5.797060| 0.0000840| 0.0092611| 0.0077715|SYT1/STX1A/RAB27A/SV2B/SYNGR1/SYP/CSPG5                                                             |     7|
## |GO:0006836 |GO:0006836 |neurotransmitter transport                      |7/47      |161/6858 |  0.0434783|       6.344126| 5.699805| 0.0000984| 0.0095836| 0.0080421|SLC17A7/SYT1/STX1A/SV2B/SYP/CSPG5/SLC1A2                                                            |     7|
## |GO:1901699 |GO:1901699 |cellular response to nitrogen compound          |10/47     |351/6858 |  0.0284900|       4.157120| 5.043880| 0.0001027| 0.0095836| 0.0080421|ACHE/ACTN2/ADCY5/LY6H/ATP2B1/ACE/GNAL/PLCB1/SLC9A1/SLC1A2                                           |    10|
## |GO:0043269 |GO:0043269 |regulation of monoatomic ion transport          |8/47      |227/6858 |  0.0352423|       5.142375| 5.272097| 0.0001312| 0.0113674| 0.0095390|ACTN2/ATP2B1/ACE/SCN4B/HOMER1/DPP10/THY1/SLC9A1                                                     |     8|
## |GO:0034762 |GO:0034762 |regulation of transmembrane transport           |8/47      |232/6858 |  0.0344828|       5.031548| 5.189203| 0.0001526| 0.0122610| 0.0102889|ACTN2/ACE/SCN4B/HOMER1/DPP10/THY1/SLC9A1/SLC1A2                                                     |     8|
## |GO:1903530 |GO:1903530 |regulation of secretion by cell                 |9/47      |300/6858 |  0.0300000|       4.377447| 4.969057| 0.0001632| 0.0122610| 0.0102889|ACHE/ADCY5/RAB27B/SYT1/STX1A/RAB27A/SV2B/PLCB1/CSPG5                                                |     9|
## |GO:1903861 |GO:1903861 |positive regulation of dendrite extension       |3/47      |17/6858  |  0.1764706|      25.749687| 8.486819| 0.0001918| 0.0122610| 0.0102889|CPNE5/SYT1/RASAL1                                                                                   |     3|
## |GO:0034329 |GO:0034329 |cell junction assembly                          |9/47      |312/6858 |  0.0288462|       4.209084| 4.819266| 0.0002194| 0.0122610| 0.0102889|ACHE/ACTN2/ACTN1/ICAM5/ACE/PLXNA4/THY1/PLXNA1/SLC9A1                                                |     9|
## |GO:0060560 |GO:0060560 |developmental growth involved in morphogenesis  |6/47      |129/6858 |  0.0465116|       6.786739| 5.511427| 0.0002249| 0.0122610| 0.0102889|OLFM1/CPNE5/SYT1/RASAL1/PLXNA4/PLXNA1                                                               |     6|
## |GO:1903859 |GO:1903859 |regulation of dendrite extension                |3/47      |18/6858  |  0.1666667|      24.319149| 8.228704| 0.0002291| 0.0122610| 0.0102889|CPNE5/SYT1/RASAL1                                                                                   |     3|

no.sig <- markers.enrich %>% 
  lapply(getElement, "result") %>% 
  sapply(\(x) sum(x$p.adjust <= 0.05))

markers.enrich %>% 
  .[no.sig > 0] %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))

## [[1]]

## 
## [[2]]

## 
## [[3]]

DISEASE-SPECIFIC PROTEINS, NO IMPUTATION

Prepare

dat.clean@colData@listData$area <- dat.clean@colData@listData$condition
dat.clean@colData@listData$areaRep <- dat.clean@colData@listData$replicate
dat.clean@colData@listData$condition <- dat.clean@colData@listData$Diagnosis
dat.clean@colData@listData$replicate <- dat.clean@colData@listData$diagRep
dat.clean@colData@rownames <- dat.clean$ID

dat.all <- dat.clean

FC

Extract results

min.l2fc = 0.1
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "FC"]

all_contrasts <- test_diff(dat.clean, type = "all")

## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'

## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL

## Warning: Partial NA coefficients for 87 probe(s)

comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))
  
res <- dep_all@elementMetadata@listData %>%
  as.data.frame()

res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
  setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj"))) %>% 
                    filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)

## $V1

## 
## $V2

## 
## $V3

## 
## $V4

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")

## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.

## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Clustering of DEPs per sample

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
plot_dist(significant = T, pal = "YlOrRd")

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
plot_heatmap(type = "centered", kmeans = F, 
             col_limit = 2, show_row_names = FALSE,
             indicate = c("condition"))

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
# 
# rowIds <- ComplexHeatmap::row_order(ht) %>% 
#   lapply(\(x) protIds[x])
# 
# go.res <- rowIds %>% 
#   lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
# 
# go.res %>% 
#   names() %>% 
#   lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))

Clustering of DEPs per comparison

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
plot_heatmap(type = "contrast", kmeans = TRUE, 
             k = 3, col_limit = 2, show_row_names = FALSE)

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

## Warning: Cannot perform kmeans clustering with missing values

Onthology on top significant proteins

markers.enrich <- res.split %>%
  lapply(dplyr::slice, seq(50)) %>% 
  .[sapply(., nrow) > 0] %>% 
  lapply(pull, name) %>%
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")

markers.enrich %<>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

markers.enrich %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))

## [[1]]

## 
## [[2]]

## 
## [[3]]

CER

Extract results

min.l2fc = 0.1
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "CER"]

all_contrasts <- test_diff(dat.clean, type = "all")

## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'

## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL

## Warning: Partial NA coefficients for 97 probe(s)

comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))
  
res <- dep_all@elementMetadata@listData %>%
  as.data.frame()

res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
  setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj"))) %>% 
                    filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)

## $V1

## 
## $V2

## 
## $V3

## 
## $V4

## Warning: ggrepel: 43 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")

## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.

## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Clustering of DEPs per sample

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
  plot_heatmap(type = "centered", kmeans = F, 
             col_limit = 2, show_row_names = FALSE,
             indicate = c("condition"))

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
# 
# rowIds <- ComplexHeatmap::row_order(ht) %>% 
#   lapply(\(x) protIds[x])
# 
# go.res <- rowIds %>% 
#   lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
# 
# go.res %>% 
#   names() %>% 
#   lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))

Clustering of DEPs per comparison

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
  plot_heatmap(type = "contrast", kmeans = F, 
             col_limit = 3, show_row_names = FALSE)

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

Onthology on top significant proteins

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## --> No gene can be mapped....

## --> Expected input gene ID: PRMT1,AURKC,MEF2A,RTEL1,TP53BP1,ANKLE1

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: NOP53,WDR48,RMI2,STOX1,CGAS,PRIMPOL

## --> return NULL...

res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## --> No gene can be mapped....

## --> Expected input gene ID: TFRC,MSH2,MSH3,SESN2,CHEK1,KPNA2

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: ENDOG,YEATS4,SPIDR,NDFIP1,NSD2,TGFB1

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: MAGEF1,RRM1,H1-4,TEX15,VPS72,TFRC

## --> return NULL...

res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))

## [[1]]

## 
## [[2]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))

## [[1]]

## 
## [[2]]

STR

Extract results

min.l2fc = 1.65
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "STR"]

all_contrasts <- test_diff(dat.clean, type = "all")

## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'

## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL

## Warning: Partial NA coefficients for 165 probe(s)

comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))
  
res <- dep_all@elementMetadata@listData %>%
  as.data.frame()

res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
  setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj"))) %>% 
                    filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)

## $V1

## Warning: ggrepel: 30 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V2

## 
## $V3

## 
## $V4

## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V5

## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4, indicate = "condition")$data %>% 
  ggplot(aes(PC1, PC2, shape = condition, col = subtype)) +
  geom_point(size = 3) +
  theme_bw()

Clustering of DEPs per sample

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
  plot_heatmap(type = "centered", kmeans = F, 
             col_limit = 3, show_row_names = F,
             indicate = c("condition"))

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

## Warning in min(x): no non-missing arguments to min; returning Inf

## Warning in max(x): no non-missing arguments to max; returning -Inf

## Warning in min(x): no non-missing arguments to min; returning Inf

## Warning in max(x): no non-missing arguments to max; returning -Inf

# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
# 
# rowIds <- ComplexHeatmap::row_order(ht) %>% 
#   lapply(\(x) protIds[x])
# 
# go.res <- rowIds %>% 
#   lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
# 
# go.res %>% 
#   names() %>% 
#   lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))

Clustering of DEPs per comparison

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
  plot_heatmap(, type = "contrast", kmeans = F, 
             col_limit = 5, show_row_names = FALSE)

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

Onthology on top significant proteins

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## --> No gene can be mapped....

## --> Expected input gene ID: KDM1A,SMCHD1,FIGNL1,OOEP,RPL24,TERF2IP

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: CHCHD4,POLG2,STOX1,RPS5,DNAJA3,RPA2

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: RIF1,H1-10,XRCC5,AKT3,NEURL4,KPNA2

## --> return NULL...

res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## --> No gene can be mapped....

## --> Expected input gene ID: TFRC,CDCA8,AURKB,H1-5,ERCC2,MRM2

## --> return NULL...

res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))

## [[1]]

## 
## [[2]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

STRINGdb network

Networks are shown for all evidence of interaction, not just physical evidence (stronger).

sdb <- STRINGdb::STRINGdb(species = 9606, score_threshold = 400, version = "12.0")

sdb.prots <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    message(x)
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_"))) %>%
      pull(name)
  }) %>% 
  setNames(names(res.split)) %>%
  .[sapply(., length) > 1]

## PSP_vs_MSA

## PSP_vs_PD

## PSP_vs_CTRL

## MSA_vs_PD

## MSA_vs_CTRL

## PD_vs_CTRL

sdb.prots %>% 
  lapply(sdb$plot_network)

## $PSP_vs_MSA
## NULL
## 
## $MSA_vs_PD
## NULL
## 
## $MSA_vs_CTRL
## NULL
## 
## $PD_vs_CTRL
## NULL

SN

Extract results

min.l2fc = 0.5
min.padj = 0.05

dat.clean <- dat.all[, dat.all$area == "SN"]

all_contrasts <- test_diff(dat.clean, type = "all")

## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'

## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL

## Warning: Partial NA coefficients for 357 probe(s)

comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))
  
res <- dep_all@elementMetadata@listData %>%
  as.data.frame()

res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
  setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj"))) %>% 
                    filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = F)

## $V1

## 
## $V2

## 
## $V3

## Warning: ggrepel: 18 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## $V4

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 100, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")

## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.

## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4, indicate = "condition")$data %>% 
  ggplot(aes(PC1, PC2, shape = condition, col = subtype)) +
  geom_point(size = 3) +
  theme_bw()

Clustering of DEPs per sample

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
plot_heatmap(type = "centered", kmeans = F, 
             col_limit = 3, show_row_names = FALSE,
             indicate = c("condition"))

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

## Warning in min(x): no non-missing arguments to min; returning Inf

## Warning in max(x): no non-missing arguments to max; returning -Inf

# protIds <- ht@ht_list$`log2 Centered intensity`@matrix %>% rownames()
# 
# rowIds <- ComplexHeatmap::row_order(ht) %>% 
#   lapply(\(x) protIds[x])
# 
# go.res <- rowIds %>% 
#   lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
# 
# go.res %>% 
#   names() %>% 
#   lapply(\(x) enrichplot::dotplot(go.res[[x]], title = x))

Clustering of DEPs per comparison

dep_all %>% 
  .[!is.na(.@elementMetadata$significant), ] %>% 
plot_heatmap(type = "contrast", kmeans = F, 
             col_limit = 3, show_row_names = FALSE)

## Warning: Missing values in '.'. Using clustering_distance = 'gower'

Onthology on top significant proteins

res.go.up <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>%
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) > 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## --> No gene can be mapped....

## --> Expected input gene ID: DNA2,WRAP53,UBQLN4,AURKB,H1-8,MGME1

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: IL10,ZNF365,EP400,TOP3A,SHLD2,IL27RA

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: HMCES,RPS27,RAD51AP1,YEATS4,RNF8,CREBBP

## --> return NULL...

res.go.down <- res.split %>% 
  names() %>% 
  lapply(\(x) {
    res.split[[x]] %>% 
      filter(!!sym(paste(x, "significant", sep = "_")),
             !!sym(paste(x, "diff", sep = "_")) < 0) %>% 
      pull(name) %>% 
      clusterProfiler::enrichGO(OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T, ont = "BP")
  }) %>% 
  setNames(names(res.split)) %>% 
  .[!sapply(., is.null)] %>% 
  .[sapply(., nrow) > 0] %>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

## --> No gene can be mapped....

## --> Expected input gene ID: METTL4,LOC102724159,UNG,H1-3,EXOSC3,TERF2IP

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: MAGEF1,H1-1,CDCA8,PARPBP,RPS5,H1-7

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: BOP1,ACTR2,HELB,SENP3,SESN2,RAD50

## --> return NULL...

res.go.up %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.up[[x]], title = x))

## [[1]]

## 
## [[2]]

## 
## [[3]]

res.go.down %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(res.go.down[[x]], title = x))

## [[1]]

## 
## [[2]]

MED

Extract results

min.l2fc = 0.1
min.padj = 0.2

dat.clean <- dat.all[, dat.all$area == "MED"]

all_contrasts <- test_diff(dat.clean, type = "all")

## Warning in test_diff(dat.clean, type = "all"): Missing values in 'dat.clean'

## Tested contrasts: PSP_vs_MSA, PSP_vs_PD, PSP_vs_CTRL, MSA_vs_PD, MSA_vs_CTRL, PD_vs_CTRL

## Warning: Partial NA coefficients for 724 probe(s)

comb <- combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_"))

dep_all <- all_contrasts

for (cc in comb) {
  dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <- p.adjust(dep_all@elementMetadata[[paste0(cc, "_p.val")]], method = "fdr")
  
   dep_all@elementMetadata[[paste0(cc, "_significant")]] <-
    (dep_all@elementMetadata[[paste0(cc, "_p.adj")]] <= min.padj) &
    (abs(dep_all@elementMetadata[[paste0(cc, "_diff")]]) >= min.l2fc)
}

dep_all@elementMetadata$significant <- dep_all@elementMetadata %>% 
  .[, grepl("_significant", names(.))] %>% 
  getElement("listData") %>% 
  as.data.frame() %>% 
  apply(1, \(x) any(x))
  
res <- dep_all@elementMetadata@listData %>%
  as.data.frame()

res.split <- comb %>% ## Måske man burde filtrere baseret på < 70 % NAs per comparison
  setNames(lapply(., \(cc) res[, c("name", colnames(res)[grep(cc, colnames(res))])] %>%
                    arrange(!!sym(paste0(cc, "_p.adj"))) %>% 
                    filter(!is.na(!!sym(paste0(cc, "_p.adj"))))), .)

Number of significant proteins

res.split %>% 
  names() %>% 
  lapply(\(cc) res.split[[cc]] %>% 
           filter(!!sym(paste0(cc, "_p.adj")) <= min.padj, abs(!!sym(paste0(cc, "_diff"))) >= min.l2fc)) %>% 
  setNames(names(res.split)) %>% 
  lapply(nrow) %>% 
  unlist() %>% 
  {data.frame(comp = names(.), value = unname(.))} %>% 
  ggplot(aes(comp, value, fill = comp)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(panel.grid.major.y = element_blank()) +
  labs(x = "", y = "No. significant proteins") +
  scale_fill_manual(values = ggsci::pal_cosmic()(7))

Volcano plots

combn(c("PSP", "MSA", "PD", "CTRL"), 2) %>% 
  as.data.frame() %>% 
  apply(2, \(x) paste(x[1], x[2], sep = "_vs_")) %>% 
  lapply(plot_volcano, dep = dep_all, label_size = 2, add_names = TRUE, adjusted = T)

## $V1

## 
## $V2

## 
## $V3

## 
## $V4

## 
## $V5

## 
## $V6

PCA

plot_pca(dep_all, x = 1, y = 2, n = 500, point_size = 4) +
  scale_shape_manual(values = rep(20, ncol(dat.clean))) +
  ggplot2::guides(shape = "none")

## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.

## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

Onthology on top significant proteins

markers.enrich <- res.split %>%
  lapply(dplyr::slice, seq(50)) %>% 
  .[sapply(., nrow) > 0] %>% 
  lapply(pull, name) %>%
  lapply(clusterProfiler::enrichGO, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL", universe = universe, readable = T)

markers.enrich %<>% 
  .[(sapply(., \(x) {
    x@result %>% 
      filter(p.adjust <= 0.05) %>% 
      nrow()})) > 0]

markers.enrich %>% 
  names() %>% 
  lapply(\(x) enrichplot::dotplot(markers.enrich[[x]], title = x))

## [[1]]

Session

Sys.time() - tt

## Time difference of 11.8523 mins

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2025.3/lib/libmkl_gf_lp64.so.2;  LAPACK version 3.12.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Berlin
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] DEP_1.26.0                  qs_0.27.3                  
##  [3] SummarizedExperiment_1.38.1 Biobase_2.68.0             
##  [5] GenomicRanges_1.60.0        GenomeInfoDb_1.44.0        
##  [7] IRanges_2.42.0              S4Vectors_0.46.0           
##  [9] BiocGenerics_0.54.0         generics_0.1.4             
## [11] MatrixGenerics_1.20.0       matrixStats_1.5.0          
## [13] pheatmap_1.0.12             biomaRt_2.64.0             
## [15] scHelper_0.0.5              ggpubr_0.6.0               
## [17] ggrepel_0.9.6               cowplot_1.1.3              
## [19] ggplot2_3.5.2               dplyr_1.1.4                
## [21] magrittr_2.0.3             
## 
## loaded via a namespace (and not attached):
##   [1] STRINGdb_2.20.0             bitops_1.0-9               
##   [3] fs_1.6.6                    ProtGenerics_1.40.0        
##   [5] enrichplot_1.28.2           httr_1.4.7                 
##   [7] RColorBrewer_1.1-3          doParallel_1.0.17          
##   [9] ggsci_3.2.0                 tools_4.5.2                
##  [11] MSnbase_2.30.1              backports_1.5.0            
##  [13] R6_2.6.1                    DT_0.33                    
##  [15] lazyeval_0.2.2              GetoptLong_1.0.5           
##  [17] withr_3.0.2                 prettyunits_1.2.0          
##  [19] gridExtra_2.3               preprocessCore_1.70.0      
##  [21] fdrtool_1.2.18              cli_3.6.5                  
##  [23] Cairo_1.6-2                 sandwich_3.1-1             
##  [25] labeling_0.4.3              conos_1.5.2                
##  [27] sass_0.4.10                 mvtnorm_1.3-3              
##  [29] readr_2.1.5                 yulab.utils_0.2.0          
##  [31] gson_0.1.0                  DOSE_4.2.0                 
##  [33] R.utils_2.13.0              dichromat_2.0-0.1          
##  [35] MetaboCoreUtils_1.16.1      plotrix_3.8-4              
##  [37] limma_3.64.0                rstudioapi_0.17.1          
##  [39] impute_1.82.0               RSQLite_2.3.11             
##  [41] gridGraphics_0.5-1          shape_1.4.6.1              
##  [43] RApiSerialize_0.1.4         gtools_3.9.5               
##  [45] car_3.1-3                   GO.db_3.21.0               
##  [47] Matrix_1.7-3                MALDIquant_1.22.3          
##  [49] imputeLCMD_2.1              abind_1.4-8                
##  [51] R.methodsS3_1.8.2           lifecycle_1.0.4            
##  [53] yaml_2.3.10                 carData_3.0-5              
##  [55] gplots_3.2.0                qvalue_2.40.0              
##  [57] SparseArray_1.8.0           BiocFileCache_2.16.0       
##  [59] Rtsne_0.17                  grid_4.5.2                 
##  [61] blob_1.2.4                  promises_1.3.2             
##  [63] crayon_1.5.3                PSMatch_1.12.0             
##  [65] shinydashboard_0.7.3        ggtangle_0.0.6             
##  [67] lattice_0.22-7              mzR_2.38.0                 
##  [69] KEGGREST_1.48.0             magick_2.8.6               
##  [71] pillar_1.10.2               knitr_1.50                 
##  [73] ComplexHeatmap_2.24.0       fgsea_1.34.0               
##  [75] rjson_0.2.23                codetools_0.2-20           
##  [77] fastmatch_1.1-6             glue_1.8.0                 
##  [79] ggfun_0.1.8                 pcaMethods_2.0.0           
##  [81] data.table_1.17.2           MultiAssayExperiment_1.34.0
##  [83] vctrs_0.6.5                 png_0.1-8                  
##  [85] treeio_1.32.0               gtable_0.3.6               
##  [87] gsubfn_0.7                  assertthat_0.2.1           
##  [89] cachem_1.1.0                xfun_0.52                  
##  [91] S4Arrays_1.8.0              mime_0.13                  
##  [93] ncdf4_1.24                  iterators_1.0.14           
##  [95] gmm_1.8                     statmod_1.5.0              
##  [97] nlme_3.1-168                ggtree_3.16.0              
##  [99] bit64_4.6.0-1               progress_1.2.3             
## [101] filelock_1.0.3              bslib_0.9.0                
## [103] affyio_1.78.0               tmvtnorm_1.6               
## [105] KernSmooth_2.23-26          colorspace_2.1-1           
## [107] DBI_1.2.3                   tidyselect_1.2.1           
## [109] chron_2.3-62                bit_4.6.0                  
## [111] compiler_4.5.2              curl_7.0.0                 
## [113] httr2_1.2.1                 xml2_1.4.0                 
## [115] DelayedArray_0.34.1         stringfish_0.16.0          
## [117] caTools_1.18.3              scales_1.4.0               
## [119] affy_1.86.0                 rappdirs_0.3.3             
## [121] stringr_1.5.1               digest_0.6.37              
## [123] rmarkdown_2.29              XVector_0.48.0             
## [125] htmltools_0.5.8.1           pkgconfig_2.0.3            
## [127] dbplyr_2.5.1                fastmap_1.2.0              
## [129] rlang_1.1.6                 GlobalOptions_0.1.2        
## [131] htmlwidgets_1.6.4           UCSC.utils_1.4.0           
## [133] shiny_1.10.0                farver_2.1.2               
## [135] jquerylib_0.1.4             zoo_1.8-14                 
## [137] jsonlite_2.0.0              BiocParallel_1.42.0        
## [139] mzID_1.46.0                 GOSemSim_2.34.0            
## [141] R.oo_1.27.1                 Formula_1.2-5              
## [143] GenomeInfoDbData_1.2.14     ggplotify_0.1.2            
## [145] patchwork_1.3.0             Rcpp_1.0.14                
## [147] proto_1.0.0                 ape_5.8-1                  
## [149] MsCoreUtils_1.20.0          sqldf_0.4-11               
## [151] vsn_3.76.0                  leidenAlg_1.1.5            
## [153] stringi_1.8.7               MASS_7.3-65                
## [155] org.Hs.eg.db_3.21.0         plyr_1.8.9                 
## [157] parallel_4.5.2              Biostrings_2.76.0          
## [159] sccore_1.0.5                splines_4.5.2              
## [161] hash_2.2.6.3                hms_1.1.3                  
## [163] Spectra_1.18.0              circlize_0.4.16            
## [165] igraph_2.1.4                QFeatures_1.18.0           
## [167] ggsignif_0.6.4              reshape2_1.4.4             
## [169] XML_3.99-0.18               evaluate_1.0.3             
## [171] RcppParallel_5.1.10         BiocManager_1.30.25        
## [173] tzdb_0.5.0                  foreach_1.5.2              
## [175] httpuv_1.6.16               tidyr_1.3.1                
## [177] purrr_1.0.4                 clue_0.3-66                
## [179] norm_1.0-11.1               BiocBaseUtils_1.10.0       
## [181] broom_1.0.8                 xtable_1.8-4               
## [183] AnnotationFilter_1.32.0     tidytree_0.4.6             
## [185] rstatix_0.7.2               later_1.4.2                
## [187] tibble_3.2.1                clusterProfiler_4.16.0     
## [189] aplot_0.2.5                 memoise_2.0.1              
## [191] AnnotationDbi_1.70.0        cluster_2.1.8.1